Outline

Original Idea was to target JAMA - thinking one of the Research Letter article type:

Research Letters are concise, focused reports of original research. These should not exceed 600 words of text and 6 references and may include up to 2 tables or figures. Online supplementary material is only allowed for brief additional and absolutely necessary methods but not for any additional results or discussion. Research Letters may have no more than 7 authors. The text should include the full name, academic degrees, and a single institutional affiliation for each author and the email address for the corresponding author. Other persons who have contributed to the study may be indicated in an Acknowledgment, with their permission, including their academic degrees, affiliation, contribution to the study, and an indication if compensation was received for their role. Letters must not duplicate other material published or submitted for publication. In general, Research Letters should be divided into the following sections: Introduction, Methods, Results, and Discussion. They should not include an abstract, but otherwise should follow all of the guidelines in Manuscript Preparation and Submission Requirements. Letters not meeting these specifications are generally not considered.

Background on California’s 2019

What is PSPS?

Potential impacts on Health

  • what we are looking at; equity as upstream drivers of who is least able to cope with PSPS
  • what we are not looking at; emergency preparedness, health impacts related to DME

Data and methods

PSPS data

  • PG&E’s records of where and when power was shut off
  • cross referenced with reports on each event from PG&E

SDoH Data

  • HPI (Indices of disadvantage)
  • CDC’s Social Vulnerability Index (maybe?)
  • specific indicators (from American Community Survey)

Assessment

  • approximating impacted population
  • statistical tests by event

Results

Characterizing Events

  • Events varied in size and location
  • smaller events occurred more frequently (tended to affect foothills - though an early event impacted Wine Country), while larger events tended to encompass both areas.
  • important differences
    • age
    • disability
    • poverty
      • food insecurity
      • insurance?
  • smaller event impacting the foothills (Paradise, Oroville) show the most pronounced disparities in SDoH
    • these areas were impacted by the most events (duration?)

Recommendations

Acting on information

  • use data on inequities of impacted communities to prioritize re-energization
  • match resources to inequity and disadvantage
  • incorporate info on inequity (related to PSPS) into emergency preparedness and response

Improve data transparency


Draft

In 2019, following the most devastating wildfire season in the state’s history, California’s largest utiility, Pacific Gas and Electric (PG&E) began using the practice of pre-emptive deenergization of utility line to avoid setting fires. During the 2019 Fire Season in California these events, refered to as Public Safety Power Shutoffs (PSPS), were used on several occasions.

When PSPS events occur residents experience a loss of power that can last beyond the


Mapping PSPS Events

Spatial data of PSPS affected areas initially received from PG&E was cross-referenced with reports on their website.

There are some reports on PG&E’s website for which we do not have spatial data about affected areas.

Spatial data from PG&E contained boundaries for 12 separate areas, which - using layer names - were cross-referenced to reports for 7 specific de-energization events (shown on the map below, click for PG&E’s estimates start/end times and number of affected customers).

event1 <- filter(pge_events, event == "June_7-9")
event2 <- filter(pge_events, event == "September_25-27")
event3 <- filter(pge_events, event == "October_5")
event4 <- filter(pge_events, event == "October_9-12")
event5 <- filter(pge_events, event == "October_23-25")
event6 <- filter(pge_events, event == "October_26-November_1")
event7 <- filter(pge_events, event == "November_20-21")

pal2 <- colorQuantile(palette = c("#2172B4","#9DCAE1","#A4DBA1","#008836"), n=4,  domain = 0:100)


maphpi <- merge(tracts_sp, { 
          
          hpi[,.(ct10, hpi2_pctile)]
  
          
        })

    leaflet() %>%
     # addProviderTiles(providers$Esri.WorldStreetMap, group = "Street Map") %>%
     #  addProviderTiles(providers$Stamen.Terrain, group = "Terrain") %>%
     #  addProviderTiles(providers$Esri.WorldImagery, group = "Sattelite") %>%
    addProviderTiles(providers$Stamen.Toner, group = "B/W") %>%  
        addPolygons(data = maphpi,
          
        color = "#444444",
        weight = 1,
        smoothFactor = 0.1,
        fillOpacity = 0.5,
        fillColor = ~ pal2(maphpi$hpi2_pctile),
        stroke = FALSE,
        popup = paste0("This location is in the top ",round(maphpi$hpi2_pctile,0),"the percentile of all tracts in the state for the healthy places index"),
        group = "HPI"
      ) %>%
    addPolygons(
      data = st_zm(event1),
      color = "#252525",
      weight = 2,
      smoothFactor = 0.1,
      fillOpacity = 0.7,
      fillColor = "#377eb8",
      stroke = TRUE,
      label = ~event1$event,
      popup = paste0("This PSPS event spanned ", round(event1$Duration_hrs,0), " hours   \n   (beginning ", as.character(event1$start_time), " and ending at ", event1$end_time,").     \n   Ultimately PG&E estimates that it affected ",format(event1$customers, big.mark = ',')," customers."),
      group = "June 7-9 (blue)"
    ) %>%
    addPolygons(
      data = st_zm(event2),
      color = "#252525",
      weight = 2,
      smoothFactor = 0.1,
      fillOpacity = 0.7,
      fillColor = "#4daf4a",
      stroke = TRUE,
      label = ~event2$event,
      popup = paste0("This PSPS event spanned ", round(event2$Duration_hrs,0), " hours   \n   (beginning ", as.character(event2$start_time), " and ending at ", event2$end_time,").     \n   Ultimately PG&E estimates that it affected ",format(event2$customers, big.mark = ',')," customers."),
      group = "Sep 25-27 (green)"
    ) %>%
   addPolygons(
      data = st_zm(event3),
      color = "#252525",
      weight = 2,
      smoothFactor = 0.1,
      fillOpacity = 0.7,
      fillColor = "#984ea3",
      stroke = TRUE,
      label = ~event3$event,
      popup = paste0("This PSPS event spanned ", round(event3$Duration_hrs,0), " hours   \n   (beginning ", as.character(event3$start_time), " and ending at ", event3$end_time,").     \n   Ultimately PG&E estimates that it affected ",format(event3$customers, big.mark = ',')," customers."),
      group = "Oct 5 (purple)"
    ) %>%
   addPolygons(
      data = st_zm(event4),
      color = "#252525",
      weight = 2,
      smoothFactor = 0.1,
      fillOpacity = 0.7,
      fillColor = "#ffff33",
      stroke = TRUE,
      label = ~event4$event,
      popup = paste0("This PSPS event spanned ", round(event4$Duration_hrs,0), " hours   \n   (beginning ", as.character(event4$start_time), " and ending at ", event4$end_time,").     \n   Ultimately PG&E estimates that it affected ",format(event4$customers, big.mark = ',')," customers."),
      group = "Oct 9-12 (yellow)"
    ) %>%
       addPolygons(
      data = st_zm(event5),
      color = "#252525",
      weight = 2,
      smoothFactor = 0.1,
      fillOpacity = 0.7,
      fillColor = "#a65628",
      stroke = TRUE,
      label = ~event5$event,
      popup = paste0("This PSPS event spanned ", round(event5$Duration_hrs,0), " hours   \n   (beginning ", as.character(event5$start_time), " and ending at ", event5$end_time,").     \n   Ultimately PG&E estimates that it affected ",format(event5$customers, big.mark = ',')," customers."),
      group = "Oct 23-25 (brown)"
    ) %>%
    addPolygons(
      data = st_zm(event6),
      color = "#252525",
      weight = 2,
      smoothFactor = 0.1,
      fillOpacity = 0.7,
      fillColor = "#e41a1c",
      stroke = TRUE,
      label = ~event6$event,
      popup = paste0("This PSPS event spanned ", round(event6$Duration_hrs,0), " hours   \n   (beginning ", as.character(event6$start_time), " and ending at ", event6$end_time,").     \n   Ultimately PG&E estimates that it affected ",format(event6$customers, big.mark = ',')," customers."),
      group = "Oct 26 - Nov 1 (red)"
    ) %>%
       addPolygons(
      data = st_zm(event7),
      color = "#252525",
      weight = 2,
      smoothFactor = 0.1,
      fillOpacity = 0.7,
      fillColor = "#ff7f00",
      stroke = TRUE,
      label = ~event7$event,
      popup = paste0("This PSPS event spanned ", round(event7$Duration_hrs,0), " hours   \n   (beginning ", as.character(event7$start_time), " and ending at ", event7$end_time,").     \n   Ultimately PG&E estimates that it affected ",format(event7$customers, big.mark = ',')," customers."),
      group = "Nov 20-21 (orange)"
    ) %>%
    addLayersControl(
    # baseGroups = c("Street Map", "Terrain", "Sattelite", "B/W"),
      overlayGroups = c("HPI","June 7-9 (blue)", "Sep 25-27 (green)", "Oct 5 (purple)","Oct 9-12 (yellow)", "Oct 23-25 (brown)" , "Oct 26 - Nov 1 (red)",  "Nov 20-21 (orange)"),
      options = layersControlOptions(collapsed = FALSE)
                )%>%  # adds a button to return to the whole state view
      addLegend("bottomleft",
                pal = pal2,
                values = na.omit(maphpi$hpi2_pctile), 
                labels = c("bottom 25%","","","healthiest 25%"),
                opacity = 1,
                title = "HPI Percentile") %>%
      setView(lat = 39.302394, lng = -122.084003, zoom = 6) %>% #defaults the view of the original map to show the entire state
      addEasyButton(easyButton(
        icon="fa-globe", title="Zoom to State",
        onClick=JS("function(btn, map){ map.setView([37.085206, -119.540085],5); }"))) %>%
      hideGroup(c("June 7-9 (blue)", "Sep 25-27 (green)", "Oct 5 (purple)","Oct 9-12 (yellow)", "Oct 23-25 (brown)" , "Oct 26 - Nov 1 (red)",  "Nov 20-21 (orange)"))
    # adds a button to return to the whole state view
      # addLegend("bottomleft",
      #           pal = pal,
      #           values = factor(c("Moderate","High","VeryHigh"), levels = c("Moderate","High","VeryHigh")),
      #           opacity = 1,
      #           title = "Fire Hazard Severity Zones"
      #           ##lables = c("X%-X% (most vulnerable)","X%-X%","X%-X%","X%-X%","X%-X% (least"%. In this tract, the "mapTemp$def," is higher than "mapTemp),))##
      # )

Statistical Testing of Specific Indicators

In order to compare population affected and not affected by specific events, we query the American Community Survey 5-year estimates for several indicators of vulnerability.

Race

% White non-Hispanic in the different areas by event.

############ DISABILITY ##############
race <- get_acs(
  geography = "tract",
  variables = c(
    # numerator = "B01001A_001", # white alone
    # numerator = "B01001B_001", # black alone
    # numerator = "B01001C_001", # American Indian Alaska alone
    # numerator = "B01001D_001", # Asian alone
    # numerator = "B01001E_001", # API alone
    # numerator = "B01001F_001", # other  alone
    # numerator = "B01001G_001", # 2 or more races
    numerator = "B01001H_001", # white alone nonhispanic
    # numerator = "B01001I_001", # Hispanic
    denominator = "B01001_001"  # total people
  ),
  state = "CA"
) %>% data.table()


race  <- race[, .(estimate = sum(estimate),
                  moe = moe_sum(estimate = estimate, moe = moe)),
              by = .(GEOID, NAME, variable)]

try1 <- merge(race, {
  readRDS("//mnt/projects/ohe/PSPS/PSPS_yn_master_tract.RDS")
  }, allow.cartesian = T)
# try1 <- na.omit(try1)

try1 <- try1[, se:= moe/1.645]

sums <- try1[,.(sum = sum(estimate, na.rm=T), 
                se_sum = (sum(se^2, na.rm=T)^0.5)), by = .(variable, PSPS_yn, event)]

totals <- sums %>% dcast.data.table(formula = PSPS_yn + event ~ variable, value.var = "sum")

sum_ses <- sums %>% dcast.data.table(formula = PSPS_yn + event ~ variable, value.var = "se_sum") %>% .[,.(PSPS_yn, event, SE_num = numerator, SE_pop = denominator)]

final <- merge(totals, sum_ses)

final <- final[, `:=` (proportion = numerator/denominator,
                       se_prop = (1/denominator)*(SE_num^2 - (((numerator/denominator)^2)*SE_pop^2))^0.5)]

final <- final[, `:=` (LL = proportion - 1.96*se_prop,
                       UL = proportion + 1.96*se_prop)]

final <- final[, event := factor(event, levels = c("June_7-9", "September_25-27","October_5","October_9-12","October_23-25","October_26-November_1","November_20-21"))] %>% setorder(event)

final <- final[, PSPS_yn := ifelse(PSPS_yn == "No","Not Affected by PSPS", "Affected by PSPS")]

hc1 <- highchart() %>% 
  hc_chart(type = "bar") %>%  
  hc_add_series(final, "point", hcaes(y = 100*proportion, x = event, group = PSPS_yn))%>% 
  hc_add_series(final, "errorbar", hcaes(low = 100*LL, high = 100*UL, x = event, group = PSPS_yn)) %>%
  hc_add_theme(hc_theme_smpl()) %>% 
  hc_xAxis(title = list(text = "PSPS Event"), 
           type = 'category'  ,
          categories =  c("June 7-9  (14hrs, 22k)", "Sep 25-27  (26hrs, 75k)","Oct 5  (18hrs, 11k) ","Oct 9-12  (90hrs, 729k)","Oct 23-25  (52hrs, 177k)","Oct 26-Nov 1  (129hrs, 941k)","Nov 20-21  (42hrs, 49k)")
          ) %>%
  hc_yAxis(title = list(text = "% white (non-Hispanic) (with 95% CI)")) %>%
  hc_title(text = "% white and non-Hispanic",
           margin = 20, align = "left") %>%
  hc_subtitle(text = "Tracts Impacted by PSPS vs Not Impacted",
              align = "left", style = list(fontWeight = "bold"))


hc1

Poverty

% Living in in the different areas by event.

############ DISABILITY ##############
poverty <- get_acs(
  geography = "tract",
  variables = c( 
    numerator = "B05010_002", # below 1.00 ratio
    numerator = "B05010_010", # income ration 1.0 - 1.99
    denominator = "B05010_001"  # total workers
  ),
  state = "CA"
) %>% data.table()


poverty  <- poverty[, .(estimate = sum(estimate),
                  moe = moe_sum(estimate = estimate, moe = moe)),
              by = .(GEOID, NAME, variable)]


try1 <- merge(poverty, {
  readRDS("//mnt/projects/ohe/PSPS/PSPS_yn_master_tract.RDS")
  }, allow.cartesian = T)
# try1 <- na.omit(try1)

try1 <- try1[, se:= moe/1.645]

sums <- try1[,.(sum = sum(estimate, na.rm=T), 
                se_sum = (sum(se^2, na.rm=T)^0.5)), by = .(variable, PSPS_yn, event)]

totals <- sums %>% dcast.data.table(formula = PSPS_yn + event ~ variable, value.var = "sum")

sum_ses <- sums %>% dcast.data.table(formula = PSPS_yn + event ~ variable, value.var = "se_sum") %>% .[,.(PSPS_yn, event, SE_num = numerator, SE_pop = denominator)]

final <- merge(totals, sum_ses)

final <- final[, `:=` (proportion = numerator/denominator,
                       se_prop = (1/denominator)*(SE_num^2 - (((numerator/denominator)^2)*SE_pop^2))^0.5)]

final <- final[, `:=` (LL = proportion - 1.96*se_prop,
                       UL = proportion + 1.96*se_prop)]

final <- final[, event := factor(event, levels = c("June_7-9", "September_25-27","October_5","October_9-12","October_23-25","October_26-November_1","November_20-21"))] %>% setorder(event)

final <- final[, PSPS_yn := ifelse(PSPS_yn == "No","Not Affected by PSPS", "Affected by PSPS")]

hc1 <- highchart() %>% 
  hc_chart(type = "bar") %>%  
  hc_add_series(final, "point", hcaes(y = 100*proportion, x = event, group = PSPS_yn))%>% 
  hc_add_series(final, "errorbar", hcaes(low = 100*LL, high = 100*UL, x = event, group = PSPS_yn)) %>%
  hc_add_theme(hc_theme_smpl()) %>% 
  hc_xAxis(title = list(text = "PSPS Event"), 
           type = 'category'  ,
          categories =  c("June 7-9  (14hrs, 22k)", "Sep 25-27  (26hrs, 75k)","Oct 5  (18hrs, 11k) ","Oct 9-12  (90hrs, 729k)","Oct 23-25  (52hrs, 177k)","Oct 26-Nov 1  (129hrs, 941k)","Nov 20-21  (42hrs, 49k)")
          ) %>%
  hc_yAxis(title = list(text = "% below 200% poverty level (with 95% CI)")) %>%
  hc_title(text = "% below 200% poverty level",
           margin = 20, align = "left") %>%
  hc_subtitle(text = "Tracts Impacted by PSPS vs Not Impacted",
              align = "left", style = list(fontWeight = "bold"))


hc1

Health Insurance

% White non-Hispanic in the different areas by event.

############ DISABILITY ##############
insurance <- get_acs(
  geography = "tract",
  variables = c(
    numerator = "B27001_005", # male under 6 no health insurance
    numerator = "B27001_008", # male 6-18 no health insurance
    numerator = "B27001_011", # male 19-25 no health insurance
    numerator = "B27001_014", # male 26-34 no health insurance
    numerator = "B27001_017", # male 35-44 no health insurance
    numerator = "B27001_020", # male 45-54 no health insurance
    numerator = "B27001_023", # male 55-64 no health insurance
    numerator = "B27001_026", # male 65-74 no health insurance
    numerator = "B27001_029", # male 75+  no health insurance
    numerator = "B27001_033", # male under 6 no health insurance
    numerator = "B27001_036", # male 6-18 no health insurance
    numerator = "B27001_039", # male 19-25 no health insurance
    numerator = "B27001_042", # male 26-34 no health insurance
    numerator = "B27001_045", # male 35-44 no health insurance
    numerator = "B27001_048", # male 45-54 no health insurance
    numerator = "B27001_051", # male 55-64 no health insurance
    numerator = "B27001_054", # male 65-74 no health insurance
    numerator = "B27001_057", # male 75+  no health insurance
    denominator = "B27001_001"  # total people
  ),
  state = "CA"
) %>% data.table()


insurance  <- insurance[, .(estimate = sum(estimate),
                  moe = moe_sum(estimate = estimate, moe = moe)),
              by = .(GEOID, NAME, variable)]

try1 <- merge(insurance, {
  readRDS("//mnt/projects/ohe/PSPS/PSPS_yn_master_tract.RDS")
  }, allow.cartesian = T)
# try1 <- na.omit(try1)

try1 <- try1[, se:= moe/1.645]

sums <- try1[,.(sum = sum(estimate, na.rm=T), 
                se_sum = (sum(se^2, na.rm=T)^0.5)), by = .(variable, PSPS_yn, event)]

totals <- sums %>% dcast.data.table(formula = PSPS_yn + event ~ variable, value.var = "sum")

sum_ses <- sums %>% dcast.data.table(formula = PSPS_yn + event ~ variable, value.var = "se_sum") %>% .[,.(PSPS_yn, event, SE_num = numerator, SE_pop = denominator)]

final <- merge(totals, sum_ses)

final <- final[, `:=` (proportion = numerator/denominator,
                       se_prop = (1/denominator)*(SE_num^2 - (((numerator/denominator)^2)*SE_pop^2))^0.5)]

final <- final[, `:=` (LL = proportion - 1.96*se_prop,
                       UL = proportion + 1.96*se_prop)]

final <- final[, event := factor(event, levels = c("June_7-9", "September_25-27","October_5","October_9-12","October_23-25","October_26-November_1","November_20-21"))] %>% setorder(event)

final <- final[, PSPS_yn := ifelse(PSPS_yn == "No","Not Affected by PSPS", "Affected by PSPS")]

hc1 <- highchart() %>% 
  hc_chart(type = "bar") %>%  
  hc_add_series(final, "point", hcaes(y = 100*proportion, x = event, group = PSPS_yn))%>% 
  hc_add_series(final, "errorbar", hcaes(low = 100*LL, high = 100*UL, x = event, group = PSPS_yn)) %>%
  hc_add_theme(hc_theme_smpl()) %>% 
  hc_xAxis(title = list(text = "PSPS Event"), 
           type = 'category'  ,
          categories =  c("June 7-9  (14hrs, 22k)", "Sep 25-27  (26hrs, 75k)","Oct 5  (18hrs, 11k) ","Oct 9-12  (90hrs, 729k)","Oct 23-25  (52hrs, 177k)","Oct 26-Nov 1  (129hrs, 941k)","Nov 20-21  (42hrs, 49k)")
          ) %>%
  hc_yAxis(title = list(text = "% without Health Insurance (with 95% CI)")) %>%
  hc_title(text = "% without Health Insurance",
           margin = 20, align = "left") %>%
  hc_subtitle(text = "Tracts Impacted by PSPS vs Not Impacted",
              align = "left", style = list(fontWeight = "bold"))


hc1

Elderly

Below is an example for the % of people over 65 in the different areas by event.

old <- tidycensus::get_acs(geography = "tract", 
              variables = c(numerator = "B01001_020", 
                            numerator = "B01001_021",
                            numerator = "B01001_022",
                            numerator = "B01001_023",
                            numerator = "B01001_024",
                            numerator = "B01001_025",
                            numerator = "B01001_044", 
                            numerator = "B01001_045",
                            numerator = "B01001_046",
                            numerator = "B01001_047",
                            numerator = "B01001_048",
                            numerator = "B01001_049",
                            denominator = "B01001_001"), 
              state = "CA") %>% data.table()

old <- old[,.(estimate = sum(estimate), 
              moe = tidycensus::moe_sum(estimate = estimate, moe = moe)), 
           by= .(GEOID, NAME, variable)]

try1 <- merge(old, {
  readRDS("//mnt/projects/ohe/PSPS/PSPS_yn_master_tract.RDS")
  }, allow.cartesian = T)
# try1 <- na.omit(try1)

try1 <- try1[, se:= moe/1.645]

sums <- try1[,.(sum = sum(estimate, na.rm=T), 
                se_sum = (sum(se^2, na.rm=T)^0.5)), by = .(variable, PSPS_yn, event)]

totals <- sums %>% dcast.data.table(formula = PSPS_yn + event ~ variable, value.var = "sum")

sum_ses <- sums %>% dcast.data.table(formula = PSPS_yn + event ~ variable, value.var = "se_sum") %>% .[,.(PSPS_yn, event, SE_num = numerator, SE_pop = denominator)]

final <- merge(totals, sum_ses)

final <- final[, `:=` (proportion = numerator/denominator,
                       se_prop = (1/denominator)*(SE_num^2 - (((numerator/denominator)^2)*SE_pop^2))^0.5)]

final <- final[, `:=` (LL = proportion - 1.96*se_prop,
                       UL = proportion + 1.96*se_prop)]

final <- final[, event := factor(event, levels = c("June_7-9", "September_25-27","October_5","October_9-12","October_23-25","October_26-November_1","November_20-21"))] %>% setorder(event)

final <- final[, PSPS_yn := ifelse(PSPS_yn == "No","Not Affected by PSPS", "Affected by PSPS")]



hc1 <- highchart() %>% 
  hc_chart(type = "bar") %>%  
  hc_add_series(final, "point", hcaes(y = 100*proportion, x = event, group = PSPS_yn))%>% 
  hc_add_series(final, "errorbar", hcaes(low = 100*LL, high = 100*UL, x = event, group = PSPS_yn)) %>%
  hc_add_theme(hc_theme_smpl()) %>% 
  hc_xAxis(title = list(text = "PSPS Event"), 
           type = 'category'  ,
          categories =  c("June 7-9  (14hrs, 22k)", "Sep 25-27  (26hrs, 75k)","Oct 5  (18hrs, 11k) ","Oct 9-12  (90hrs, 729k)","Oct 23-25  (52hrs, 177k)","Oct 26-Nov 1  (129hrs, 941k)","Nov 20-21  (42hrs, 49k)")
          ) %>%
  hc_yAxis(title = list(text = "% over 65 (with 95% CI)")) %>%
  hc_title(text = "% of population over 65",
           margin = 20, align = "left") %>%
  hc_subtitle(text = "Tracts Impacted by PSPS vs Not Impacted",
              align = "left", style = list(fontWeight = "bold"))


hc1

Youth

% of people under age 5 in the different areas by event.

young <- tidycensus::get_acs(geography = "tract", 
              variables = c(numerator = "B01001_003", # male under 5
                            numerator = "B01001_027", # female under 5
                            denominator = "B01001_001"), 
              state = "CA") %>% data.table()

young <- young[,.(estimate = sum(estimate), 
              moe = tidycensus::moe_sum(estimate = estimate, moe = moe)), 
           by= .(GEOID, NAME, variable)]

try1 <- merge(young, {
  readRDS("//mnt/projects/ohe/PSPS/PSPS_yn_master_tract.RDS")
  }, allow.cartesian = T)
# try1 <- na.omit(try1)

try1 <- try1[, se:= moe/1.645]

sums <- try1[,.(sum = sum(estimate, na.rm=T), 
                se_sum = (sum(se^2, na.rm=T)^0.5)), by = .(variable, PSPS_yn, event)]

totals <- sums %>% dcast.data.table(formula = PSPS_yn + event ~ variable, value.var = "sum")

sum_ses <- sums %>% dcast.data.table(formula = PSPS_yn + event ~ variable, value.var = "se_sum") %>% .[,.(PSPS_yn, event, SE_num = numerator, SE_pop = denominator)]

final <- merge(totals, sum_ses)

final <- final[, `:=` (proportion = numerator/denominator,
                       se_prop = (1/denominator)*(SE_num^2 - (((numerator/denominator)^2)*SE_pop^2))^0.5)]

final <- final[, `:=` (LL = proportion - 1.96*se_prop,
                       UL = proportion + 1.96*se_prop)]

final <- final[, event := factor(event, levels = c("June_7-9", "September_25-27","October_5","October_9-12","October_23-25","October_26-November_1","November_20-21"))] %>% setorder(event)

final <- final[, PSPS_yn := ifelse(PSPS_yn == "No","Not Affected by PSPS", "Affected by PSPS")]



hc1 <- highchart() %>% 
  hc_chart(type = "bar") %>%  
  hc_add_series(final, "point", hcaes(y = 100*proportion, x = event, group = PSPS_yn))%>% 
  hc_add_series(final, "errorbar", hcaes(low = 100*LL, high = 100*UL, x = event, group = PSPS_yn)) %>%
  hc_add_theme(hc_theme_smpl()) %>% 
  hc_xAxis(title = list(text = "PSPS Event"), 
           type = 'category'  ,
          categories =  c("June 7-9  (14hrs, 22k)", "Sep 25-27  (26hrs, 75k)","Oct 5  (18hrs, 11k) ","Oct 9-12  (90hrs, 729k)","Oct 23-25  (52hrs, 177k)","Oct 26-Nov 1  (129hrs, 941k)","Nov 20-21  (42hrs, 49k)")
          ) %>%
  hc_yAxis(title = list(text = "% under 5 (with 95% CI)")) %>%
  hc_title(text = "% of population under 5",
           margin = 20, align = "left") %>%
  hc_subtitle(text = "Tracts Impacted by PSPS vs Not Impacted",
              align = "left", style = list(fontWeight = "bold"))


hc1

Disabled

% of people with a disability in the different areas by event.

############ DISABILITY ##############


disability <- tidycensus::get_acs(geography = "tract", 
              variables = c(numerator = "C18130_010", # with a disability 18-64  
                            numerator = "C18130_017",  # with a disability 65+  
                            denominator = "C18130_001"), 
              state = "CA") %>% data.table()

disability <- disability[,.(estimate = sum(estimate), 
              moe = tidycensus::moe_sum(estimate = estimate, moe = moe)), 
           by= .(GEOID, NAME, variable)]

try1 <- merge(disability, {
  readRDS("//mnt/projects/ohe/PSPS/PSPS_yn_master_tract.RDS")
  }, allow.cartesian = T)
# try1 <- na.omit(try1)

try1 <- try1[, se:= moe/1.645]

sums <- try1[,.(sum = sum(estimate, na.rm=T), 
                se_sum = (sum(se^2, na.rm=T)^0.5)), by = .(variable, PSPS_yn, event)]

totals <- sums %>% dcast.data.table(formula = PSPS_yn + event ~ variable, value.var = "sum")

sum_ses <- sums %>% dcast.data.table(formula = PSPS_yn + event ~ variable, value.var = "se_sum") %>% .[,.(PSPS_yn, event, SE_num = numerator, SE_pop = denominator)]

final <- merge(totals, sum_ses)

final <- final[, `:=` (proportion = numerator/denominator,
                       se_prop = (1/denominator)*(SE_num^2 - (((numerator/denominator)^2)*SE_pop^2))^0.5)]

final <- final[, `:=` (LL = proportion - 1.96*se_prop,
                       UL = proportion + 1.96*se_prop)]

final <- final[, event := factor(event, levels = c("June_7-9", "September_25-27","October_5","October_9-12","October_23-25","October_26-November_1","November_20-21"))] %>% setorder(event)

final <- final[, PSPS_yn := ifelse(PSPS_yn == "No","Not Affected by PSPS", "Affected by PSPS")]



hc1 <- highchart() %>% 
  hc_chart(type = "bar") %>%  
  hc_add_series(final, "point", hcaes(y = 100*proportion, x = event, group = PSPS_yn))%>% 
  hc_add_series(final, "errorbar", hcaes(low = 100*LL, high = 100*UL, x = event, group = PSPS_yn)) %>%
  hc_add_theme(hc_theme_smpl()) %>% 
  hc_xAxis(title = list(text = "PSPS Event"), 
           type = 'category'  ,
          categories =  c("June 7-9  (14hrs, 22k)", "Sep 25-27  (26hrs, 75k)","Oct 5  (18hrs, 11k) ","Oct 9-12  (90hrs, 729k)","Oct 23-25  (52hrs, 177k)","Oct 26-Nov 1  (129hrs, 941k)","Nov 20-21  (42hrs, 49k)")
          ) %>%
  hc_yAxis(title = list(text = "% with a disability (with 95% CI)")) %>%
  hc_title(text = "% of population living with a disability",
           margin = 20, align = "left") %>%
  hc_subtitle(text = "Tracts Impacted by PSPS vs Not Impacted",
              align = "left", style = list(fontWeight = "bold"))


hc1

Food Insecurity

% of Households receiving Food Stamps/SNAP in the past 12 months in the different areas by event.

############ DISABILITY ##############
snap <- get_acs(
  geography = "tract",
  variables = c(
    numerator = "B22001_002", # receiving SNAP  
    # numerator = "B22001_003", # receiving SNAP over 60
    denominator = "B22001_001"  # total
  ),
  state = "CA"
) %>% data.table()

snap  <- snap[, .(estimate = sum(estimate),
                               moe = moe_sum(estimate = estimate, moe = moe)),
                           by = .(GEOID, NAME, variable)]

try1 <- merge(snap, {
  readRDS("//mnt/projects/ohe/PSPS/PSPS_yn_master_tract.RDS")
  }, allow.cartesian = T)
# try1 <- na.omit(try1)

try1 <- try1[, se:= moe/1.645]

sums <- try1[,.(sum = sum(estimate, na.rm=T), 
                se_sum = (sum(se^2, na.rm=T)^0.5)), by = .(variable, PSPS_yn, event)]

totals <- sums %>% dcast.data.table(formula = PSPS_yn + event ~ variable, value.var = "sum")

sum_ses <- sums %>% dcast.data.table(formula = PSPS_yn + event ~ variable, value.var = "se_sum") %>% .[,.(PSPS_yn, event, SE_num = numerator, SE_pop = denominator)]

final <- merge(totals, sum_ses)

final <- final[, `:=` (proportion = numerator/denominator,
                       se_prop = (1/denominator)*(SE_num^2 - (((numerator/denominator)^2)*SE_pop^2))^0.5)]

final <- final[, `:=` (LL = proportion - 1.96*se_prop,
                       UL = proportion + 1.96*se_prop)]

final <- final[, event := factor(event, levels = c("June_7-9", "September_25-27","October_5","October_9-12","October_23-25","October_26-November_1","November_20-21"))] %>% setorder(event)

final <- final[, PSPS_yn := ifelse(PSPS_yn == "No","Not Affected by PSPS", "Affected by PSPS")]



hc1 <- highchart() %>% 
  hc_chart(type = "bar") %>%  
  hc_add_series(final, "point", hcaes(y = 100*proportion, x = event, group = PSPS_yn))%>% 
  hc_add_series(final, "errorbar", hcaes(low = 100*LL, high = 100*UL, x = event, group = PSPS_yn)) %>%
  hc_add_theme(hc_theme_smpl()) %>% 
  hc_xAxis(title = list(text = "PSPS Event"), 
           type = 'category'  ,
          categories =  c("June 7-9  (14hrs, 22k)", "Sep 25-27  (26hrs, 75k)","Oct 5  (18hrs, 11k) ","Oct 9-12  (90hrs, 729k)","Oct 23-25  (52hrs, 177k)","Oct 26-Nov 1  (129hrs, 941k)","Nov 20-21  (42hrs, 49k)")
          ) %>%
  hc_yAxis(title = list(text = "% receiving food stamps (with 95% CI)")) %>%
  hc_title(text = "% of Households receiving Food Stamps (past 12 months)",
           margin = 20, align = "left") %>%
  hc_subtitle(text = "Tracts Impacted by PSPS vs Not Impacted",
              align = "left", style = list(fontWeight = "bold"))


hc1

Vehicles

% without vehicle available in the different areas by event.

############ DISABILITY ##############
vehicles <- get_acs(
  geography = "tract",
  variables = c(
    # numerator = "B22001_002", # receiving SNAP  
    numerator = "B08141_002", # workers without a vehicle
    denominator = "B08141_001"  # total workers
  ),
  state = "CA"
) %>% data.table()


vehicles  <- vehicles[, .(estimate = sum(estimate),
                  moe = moe_sum(estimate = estimate, moe = moe)),
              by = .(GEOID, NAME, variable)]


try1 <- merge(vehicles, {
  readRDS("//mnt/projects/ohe/PSPS/PSPS_yn_master_tract.RDS")
  }, allow.cartesian = T)
# try1 <- na.omit(try1)

try1 <- try1[, se:= moe/1.645]

sums <- try1[,.(sum = sum(estimate, na.rm=T), 
                se_sum = (sum(se^2, na.rm=T)^0.5)), by = .(variable, PSPS_yn, event)]

totals <- sums %>% dcast.data.table(formula = PSPS_yn + event ~ variable, value.var = "sum")

sum_ses <- sums %>% dcast.data.table(formula = PSPS_yn + event ~ variable, value.var = "se_sum") %>% .[,.(PSPS_yn, event, SE_num = numerator, SE_pop = denominator)]

final <- merge(totals, sum_ses)

final <- final[, `:=` (proportion = numerator/denominator,
                       se_prop = (1/denominator)*(SE_num^2 - (((numerator/denominator)^2)*SE_pop^2))^0.5)]

final <- final[, `:=` (LL = proportion - 1.96*se_prop,
                       UL = proportion + 1.96*se_prop)]

final <- final[, event := factor(event, levels = c("June_7-9", "September_25-27","October_5","October_9-12","October_23-25","October_26-November_1","November_20-21"))] %>% setorder(event)

final <- final[, PSPS_yn := ifelse(PSPS_yn == "No","Not Affected by PSPS", "Affected by PSPS")]

hc1 <- highchart() %>% 
  hc_chart(type = "bar") %>%  
  hc_add_series(final, "point", hcaes(y = 100*proportion, x = event, group = PSPS_yn))%>% 
  hc_add_series(final, "errorbar", hcaes(low = 100*LL, high = 100*UL, x = event, group = PSPS_yn)) %>%
  hc_add_theme(hc_theme_smpl()) %>% 
  hc_xAxis(title = list(text = "PSPS Event"), 
           type = 'category'  ,
          categories =  c("June 7-9  (14hrs, 22k)", "Sep 25-27  (26hrs, 75k)","Oct 5  (18hrs, 11k) ","Oct 9-12  (90hrs, 729k)","Oct 23-25  (52hrs, 177k)","Oct 26-Nov 1  (129hrs, 941k)","Nov 20-21  (42hrs, 49k)")
          ) %>%
  hc_yAxis(title = list(text = "% without a vehicle (with 95% CI)")) %>%
  hc_title(text = "% of without access to a vehicle",
           margin = 20, align = "left") %>%
  hc_subtitle(text = "Tracts Impacted by PSPS vs Not Impacted",
              align = "left", style = list(fontWeight = "bold"))


hc1